Example Usage

import chalc as ch
import numpy as np
import matplotlib.pyplot as plt

For our data we sample 100 points on a circle with some noise and 100 points from inside the unit disk.

rng = np.random.default_rng(40)
num_points_circle = 100
num_points_disk = 100
mean = [0, 0]
cov = np.eye(2)*0.01
circle = np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_circle)]).T +\
    rng.multivariate_normal(mean, cov, num_points_circle).T # points as columns
disk = np.sqrt(rng.random(num_points_disk)) * np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_disk)]).T
points = np.concatenate((circle, disk), axis=1)
plt.scatter(circle[0, :], circle[1, :], s=10)
plt.scatter(disk[0, :], disk[1, :], s=10)
plt.gca().set_aspect('equal')
colours = [0]*num_points_circle + [1]*num_points_disk
plt.title('Point cloud')
plt.show()

We can compute the chromatic alpha/Delaunay–Cech/Delaunay–Rips complex $K$ of the point cloud using the functions chromatic.alpha, chromatic.delcech, and chromatic.delrips respectively. In each of these cases, \(K\) has far fewer simplices than either the Čech or Vietoris–Rips complex, which each have \(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices.

K, _ = ch.chromatic.alpha(points, colours)
printmd(f'$K$ has {len(K.simplices[1])} 1-simplices and {len(K.simplices[2])} 2-simplices')

\(K\) has 954 1-simplices and 1312 2-simplices

\(K\) is a FilteredComplex object, and simplices in the filtration are objects of type Simplex (although you should not generally need to interact with the API of these objects). We can visualise the filtration at any given time using plotting.draw_filtration. For example, at time 0.3 we have:

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
ch.plotting.draw_filtration(K, points, time=0.3, include_colours=[0], ax=ax[0])
ax[0].set_title('Colour 0')
ch.plotting.draw_filtration(K, points, time=0.3, include_colours=[1], ax=ax[1])
ax[1].set_title('Colour 1')
ch.plotting.draw_filtration(K, points, time=0.3, ax=ax[2])
ax[2].set_title('Both colours')
for i in range(3):
    ax[i].set_xlim(-1.2, 1.2)
    ax[i].set_ylim(-1.2, 1.2)
fig.suptitle('Complexes at $t=0.3$')
plt.tight_layout()
plt.show()

The functions sixpack.from_filtration and sixpack.compute compute the 6-pack of persistent homology diagrams of the chromatic alpha/Delaunay–Čech complex/Delaunay–Rips complexes.

# using the previously computed filtration
dgms_alpha = ch.sixpack.from_filtration(K, dom=[0])
# or directly from the point cloud
dgms_delcech = ch.sixpack.compute(points, colours, dom=[0], method="delcech")
dgms_delrips = ch.sixpack.compute(points, colours, dom=[0], method="delrips")

Each 6-pack is an object of type DiagramEnsemble, with attributes ker, cok, im, dom, cod, rel, entrance_times, and dimensions. Each of the first 6 of those attributes is an object of type SimplexPairings that has sets of paired and unpaired simplices, while the last two attributes are lists indicate the filtration time and dimension of those simplices. For example-

print(dgms_alpha.ker)

Paired: [(1187, 1200), (1445, 1453), (358, 409), (436, 658), (584, 883), (1697, 1701), (1995, 1999), (228, 235), (333, 581), (361, 456), (454, 471), (1648, 1650), (1486, 1499), (1447, 1495), (1461, 1494), (1565, 1600), (429, 502), (1510, 1514), (606, 626), (1248, 1252), (1267, 1268), (387, 587), (487, 531), (506, 596), (1454, 1457), (292, 383), (348, 528), (1329, 1331), (421, 491), (252, 438), (1592, 1671), (280, 317), (1366, 1408), (359, 412), (1464, 1470), (2013, 2433), (345, 376), (282, 338), (413, 428), (544, 646), (352, 540), (262, 372), (1692, 1693), (1360, 1458), (1393, 1403), (1505, 1516), (331, 416)]
Unpaired: []

A convenience function get is provided to extract the diagrams as a matrix of birth/death times:

np.set_printoptions(threshold=10)
print(dgms_alpha.get('ker', 1)) # get the kernel in dimension one
# dgms_alpha.get('ker', [0, 1]) to get the kernel in dimensions zero and one
# dgms_alpha.get('ker') to get the kernel in all dimensions from zero to max(dgms_alpha.dimensions)

[[0.0566658  0.06077279]
 [0.10478006 0.10515716]
 [0.14259426 0.14317964]
 ...
 [0.0909663  0.10600529]
 [0.09605532 0.09782018]
 [0.1135806  0.1155313 ]]

The convenience functions plotting.plot_sixpack and plotting.plot_diagram are provided for plotting either a 6-pack or a specific diagram from a 6-pack.

fig1, ax1 = ch.plotting.plot_sixpack(dgms_alpha, tolerance = 1e-3)
fig1.suptitle('Chromatic Alpha')
fig1.set_figwidth(15)
fig1.set_figheight(9)
fig1.subplots_adjust(top=0.92)
plt.show()
fig2, ax2 = plt.subplots(1, 3)
# You can specify the dimensions of features to include in the plot.
ch.plotting.plot_diagram(dgms_alpha, 'rel', dimensions = {0, 2}, truncation = 0.3, ax = ax2[0], tolerance = 1e-3)
ax2[0].set_title('Chromatic Alpha')
# You can also specify a single dimension as an integer.
ch.plotting.plot_diagram(dgms_delcech, 'rel', dimensions = 1, truncation = 0.3, ax = ax2[1], tolerance = 1e-3)
ax2[1].set_title('Chromatic Delaunay-Cech')
# If no dimensions are specified, all features will be included.
ch.plotting.plot_diagram(dgms_delrips, 'rel', truncation = 0.3, ax = ax2[2], tolerance = 1e-3)
ax2[2].set_title('Chromatic Delaunay-Rips')
fig2.suptitle('Individual relative diagrams')
fig2.set_figwidth(15)
plt.show()

We can save the diagrams to file or load a diagram from file:

with h5py.File('test.h5', 'w') as f:
    ch.sixpack.save_diagrams(dgms_alpha, f)

with h5py.File('test.h5', 'r') as f:
    dgms_alpha = ch.sixpack.load_diagrams(f)

We can visualise the 2-skeleton of the filtration for points in 2D:

animation = ch.plotting.animate_filtration(
    K, points, filtration_times=list(np.linspace(0, 1.0, 45)), animation_length=5)
animation

Once Loop Reflect